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At the mesoscale friction occurs through the breaking and formation of local contacts. This is 
often described by the earthquake-like model which requires numerical studies. We show that this 
phenomenon can also be described by a master equation, which can be solved analytically in some 
cases and provides an efficient numerical solution for more general cases. We examine the effect 
of temperature and aging of the contacts and discuss the statistical properties of the contacts for 
different situations of friction and their implications, particularly regarding the existence of stick- 
slip. 
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I. INTRODUCTION 

In spite of its crucial practical importance, friction is 
still far from being fully explained [iHS] ■ Besides a proper 
explanation of the static friction laws, the dynamical as- 
pects of friction are even less understood. This is exem- 
plified by the problem posed by the lack of a well estab- 
lished mechanism for the familiar stick-slip phenomenon 
that one perceives with a door's creak or the playing 
of a violin with a bow. Many phenomena in nature, 
where one part of a system moves in contact with an- 
other part, exhibit such a stick-slip motion which changes 
to smooth sliding with the increase of the driving ve- 
locity Microscopically this phenomenon was first 
explained by Robbins and Thompson [5] by a melting- 
freezing mechanism: a thin lubricant film between the 
moving surfaces melts during slips and solidifies again 
at sticks. Such a behavior is typical for conventional, 
or "soft" lubricants, when the lubricant-surface interac- 
tion is stronger than the interaction between the lubri- 
cant molecules |3J. A "hard" lubricant, which remains 
in a solid state during slips, also demonstrates the stick- 
slip to smooth sliding transition which now emerges due 
to inertia [!, @. In both cases, however, the transi- 
tion from stick-slip to smooth sliding in a planar tribo- 
logical contact is found to occur at a driving velocity 
v c ~ 10~ 2 c ~ 1 — 10 m/s (where c is the sound velocity), 
which is more than six orders of magnitude higher than 
experimentally observed [1, 0, @] ■ This leads to the con- 
clusion that microscopic mechanisms of stick-slip have 
little relevance at the macroscopic level. 

At the macroscopic scale, the stick-slip and smooth- 
sliding regimes can be explained with a phenomenolog- 
ical theory d H i, G3- 

It is based on a "contact-age 
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function" which depends on the previous history of the 
system. This theory leads to an excellent agreement with 
experiments if the model parameters are suitably chosen. 
Unfortunately, this approach remains purely phenomeno- 
logical. The corresponding equations cannot be derived 
from a microscopic-scale analysis. 

Another explanation is based on earthquake-like mod- 
els (the EQ model) such as the Burridge-Knopoff spring- 
block model [ll|, U3l later developed by Olami, Feder, 
and Christensen Il3l . It was adjusted to describe fric- 
tion by Persson [14j | and then explored in a number of 
works BUM!- In this model the two solid blocks touch 
one another at point junctions, corresponding to asper- 
ities of the surfaces which pin the relative position of 
the blocks (Fig. [IJ. A local force, exerted by the mov- 
ing top block, is associated to each contact, as well as 
a threshold for breaking. The model describes friction 
as resulting from the combined effects of the breaking 
of contacts and their pinning again after the release of 
their internal stress. Computer simulations of a simpli- 
fied variant of the EQ model, where all contacts have the 
same threshold 0, Hj] , showed that the earthquake-like 
model reproduces the experimentally observed transition 
from stick-slip to smooth sliding, provided the follow- 
ing assumptions are made: (?) the model must be two- 
dimensional, (ii) the spatial distribution of contacts must 
be random, (Hi) it should exist an interaction between 
the contacts, and (iv) the model must incorporate an 
increase of the static threshold on the time of station- 
ary contact similarly to the "age-function" assumption 
of the phenomenological theory. Recently 0, [l?} it be- 
came clear that a crucial ingredient of the EQ-like model 
in tribological applications is the distribution of static 
yield thresholds for the breaking of individual contacts. 
In the early variants of the model, all contacts had been 
assumed to be identical for the sake of simplicity, and 
a distribution of thresholds appeared only implicitly due 
to temperature fluctuations [14j or due to interaction be- 
tween the contacts As a matter of fact, the EQ model 
with identical contacts is a singular case [l6|. It admits 
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a periodic solution which can be interpreted as a form 
of stick-slip. Nonetheless this solution remains largely 
unphysical since for example it ceases to exist as soon 
as nonequivalent contacts are considered, whatever their 
precise properties. As soon as a finite width distribution 
of yield thresholds is taken into account, the solution of 
the model in the quasi-static limit always approaches a 
physical solution with smooth sliding [161 117] . Incorpo- 
rating a threshold distribution that does not evolve by 
the breaking and re-formation of the contacts allows one 
to find the steady-state solution of the EQ model ana- 
lytically, and more importantly to find conditions for ap- 
pearance of the elastic instability, which is the necessary 
condition for the stick-slip to emerge [!, [13] • 

The approach based on the earthquake-like model has, 
however, two weak points. First, most results can only 
be obtained with computer simulation. Second, the na- 
ture of contacts is not clearly specified. While for two 
rough solid surfaces the contacts may be associated with 
real asperities, for two ideal flat mica surfaces in the sur- 
face force apparatus (SFA) experiments [2C§ the nature 
of contacts is not clear. Persson [14{ supposed that they 
may correspond to "solid islands" in the fluidized lubri- 
cant, but such an assumption remains on a speculative 
level. 

Therefore, whether one looks at friction from the 
macroscopic or from the microscopic viewpoint, the theo- 
retical understanding is not satisfactory. What is lacking 
is a theory that would not separately consider the two 
extreme limits, macroscopic or microscopic, but couple 
them in a mesoscale approach. Our goal in the present 
work is to explore an intermediate approach, which starts 
from the properties of individual contacts and deduces 
macroscopic laws from an analytical description based on 
a statistical analysis. It is interesting to notice that the 
introduction of statistics in the theory of friction [2l| al- 
ready allowed a significant progress in the understanding 
of static friction, but here we are concerned with dynam- 
ical aspects which enter through the continuous breaking 
and re-forming of many local contacts. We introduce a 
master equation (ME) which describes this phenomenon 
and couples local events with macroscopic properties. It 
can be solved analytically in cases which are particularly 
relevant, or studied numerically in more complex situ- 
ations, much more efficiently than with simulations of 
discrete EQ models. Its major interest is to split the 
analysis in two independent parts: (i) the calculation of 
the friction force given by the master equation provided 
the statistical properties of the contacts are known, and 
(it) the study of the statistical properties of the contacts, 
which may need inputs from the microscopic scale. 

A preliminary report of the results have been presented 
in a short letter [17| • In this paper we discuss the mas- 
ter equation approach more thoroughly and explore some 
of the consequences of this new viewpoint on friction 
by examining several issues such as temperature effects 
fSec. HV|) and the aging of the contacts (Sec.|V]). We dis- 
cuss the origin of the statistical properties of the contacts 



(Sec. IVI[) . The paper begins by a brief review of some 
results of the EQ model (Sec. [IT]) and the presentation of 
the master equation approach and some of its analytical 
solutions in Sec. IIIII 



II. EARTHQUAKE-LIKE MODEL 

Let us begin with the earthquake-like model, which be- 
longs to the class of cellular automaton models. This sec- 
tion, which follows earlier works [14| , provides a reference 
for the master equation approach. We use a variant of 
the Burridge-Knopoff spring-block model of earthquakes 
adapted to tribology problems by Persson [14j . The con- 
tacts form an array and a local force fi(t) is associated 
with each contact. All contacts are connected through 
springs of strength ki , corresponding to their shear elastic 
constant, with the top block moving with a velocity v and 
coupled frictionally with the fixed bottom block as shown 
in Fig. [TJ The contacts may also interact elastically be- 
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FIG. 1: (Color online): The earthquake-like model. Friction 
occurs between the lower plate and the top block schematized 
by two plates connected by a network of springs to describe 
its internal shear stress. Thr contacts are represented by the 
disks on the lower plate and their elasticity is modeled by 
the springs connecting the lower plate and the bottom of the 
top block. The arrow indicate the displacement of the block. 
X and Xt are the coordinates of the bottom and top of the 
block. 

tween themselves. As the top block moves, the surface 
stress at any contact increases, fi(t) = kiXi{t), where 
Xi(t) is the shift of the ith junction from its non-stressed 
position. A single contact is assumed to be pinned whilst 
fi{t) < fsi, where f 8 i is the static friction threshold for 
the given contact. When the force reaches f S i, a rapid 
local slip takes place, during which the local stress in 
the block drops to the value ju ~ while the elastic 
energy stored in the junction is released in the form of 
phonons into the bulk. Then the contact is pinned again, 
and the whole process repeats itself. Let N c be the num- 
ber of contacts (asperities, bridges, etc.) at the inter- 
face. Each contact is characterized by its area Ai. The 
value of fcj, according to Persson [14|, can be estimated 
as ki ~ pc 2 y/A~i, where p is the mass density and c is the 
transverse sound velocity of the material which forms the 
contacts (see Appendix [AJ . Let us assume that the sub- 
strates are rigid, i.e., the elastic constant of the substrate 
is infinite, K = oo, and also neglect elastic interactions 
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between the contacts through the substrate as it was done 
in the statistical studies of static friction [21| . 

Let P c (x s ) be the normalized probability distribution 
of values of the thresholds x S i = f s i/ki at which contacts 
break, where f s i is the maximum force that contact i can 
sustain. If we denote by x s its average value and by a s 
its standard deviation, a typical example is the Gaussian 
P c (x) = G(x;x s ,a s ), where 



G(x; x, a)) 



1 



exp 



a 



(1) 



In numerics, the distribution P c (x) is defined on 
the interval < x < x m , where x m ^> x s 
so that we use the corrected distribution P c (x) — 
N{P G {x) - P G (0) - {x/x m ) [P G {x m ) - P G (0)]}, where 
XT is the normalization constant, so that P c (x) satisfies 
the condition P c (0) = P c (x m ) = 0. 

To describe the kinetics of the model, we introduce the 
distribution Q(x;X) of the stretchings Xi when the top 
substrate is at a position X. It is normalized by 



dxQ(x;X) = 1. 



(2) 



Let all contacts start from an initial distribution 
Q(x;0) = Qmi(x) corresponding to the Gaussian one, 
Q ini (x) = G(x;x in i,CTi n i) with x ini < x s and cr in i < a s . 
Let us now apply an adiabatically increasing force F to 
the top substrate, while the bottom substrate remains 
fixed. The force will induce a displacement X of the 
top substrate. According to third Newton's law (the law 
of action and reaction), in the adiabatic regime of quasi- 
equilibrium F must be compensated by the sum of elastic 
forces in the contacts, 



F(X) = N c (ki) dxxQ(x;X) 



(3) 



where we assumed that ki and Xi are independent random 
variables. 

With the increase of F and consequently X, the 
stretching of a contact i grows by the value X with re- 
spect to its initial value, until it reaches the threshold 
stretching x S i for the given contact. At this point the 
contact slides and Xi drops; we assume that the sliding 
is rapid and that the stretching drops to zero, Xbi = 0. 
In the numerical algorithm, the contact closest to the 
threshold is found and the displacement necessary to pro- 
voke a slip in this contact is added to all contacts. When 
the contact i reaches the threshold and Xi is set to zero, a 
new threshold is assigned to this contact randomly from 
the distribution P c {x). Then the whole process repeats 
itself. 

The evolution of the system is shown in Fig. [2] for short 
times and in Fig. [3] for long times. One can see that, as 
time goes, the initial distribution approaches a station- 
ary distribution Q s (x), where the total force F becomes 
independent on X, and the final distribution (see Fig. 0} 
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FIG. 2: (a) Short-time evolution of the EQ model. Solid 
curves show the distribution Q(x; X) for incrementally in- 
creasing values of X (with the step AX = 0.05). The distri- 
bution P c (x) is Gaussian with x s = 1 and a s — 0.05, the 
initial distribution Qi n i(x) is Gaussian with Xini = and 
Cini = 0.025. For clarity each graph of Q(x;X) is shifted 
by 0.025 along the x direction and by —1 along the Q direc- 
tion with respect to the previous one. (b) The corresponding 
dependence F(X). The elastic constants are (ki) — 1 and the 
top block is assumed to be fully rigid. 



is independent of the initial one. An elegant mathemati- 
cal proof of this statement for a simplified version of the 
EQ model was presented in Ref. [Tg] ; the statement is 
valid for any distribution P c (x) except the singular case 
of P c {x) = 5{x — x s ). 

Note that the probability distribution P c (x), which de- 
termines the value of the static threshold x S i for new- 
born contacts is different from the concrete realization 
of the distribution of static thresholds P xs (x) (i.e., the 
histogram calculated over the array {x S i}). While at 
the beginning P xs {x) = P c {x), then the function P xs (x) 
changes with time and finally takes the form shown in 
Fig. 2] (lower panel) by solid curve and circles. This is 
clearly demonstrated in Fig. [20] of Appendix IB1 where we 
consider a simple case of rectangularly-shaped function 

Pc{x). 



III. MASTER EQUATION APPROACH 

A. Master equation 

Let us now introduce an analytical description of the 
evolution of the distribution Q(x; X). When the top sub- 
strate is displaced by a small amount AX > (the case 
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FIG. 3: The same as in Fig. [5] for long times with the in- 
crement AX ~ 1.05. While the initial distribution extended 
to negative x values, according to the model this is not the 
case for further ones, which gives rise to what appears as a 
bump at x — (taking into account the shifts introduced for 
clarity) . 



of AX < is considered in Appendix [Cj, this increases 
all displacements of the contacts by AX too. This dis- 
placement leads to three kinds of changes in the distribu- 
tion Q(x; X): first, there is a global shift for all contacts, 
second, some contacts break because their stretchings be- 
come too large, and third, those broken contacts will form 
again, at a lower stretching after a slip at the scale of 
those contacts, which locally reduces the tension within 
the top substrate (or the tension of the corresponding 
springs in the spring-block model). These three contri- 
butions can be written as a master equation for Q(x; X): 

Q(x;X+AX) = Q(x-AX;X)-AQ-(x;X)+AQ+(x;X). 

(4) 

The first term in the r.h.s. of Eq. (fj| is just the shift. 

The second term AQ_(x; X) designates the variation 
of the distribution due to the breaking of some contacts. 
We can formally write it as 



AQ_ (a:; X) = P{x) AX Q(x; X) 



(5) 



where P(x) AX is the fraction of contacts displaced by 
x from their fully relaxed state which break when their 
displacement is increased by AX. This fraction is related 
to the individual properties of the contacts, which are de- 
termined by the distribution of static thresholds of each 
contact P c (x) defined in Sec. HU According to the def- 
inition of P c (x), the total number of unbroken contacts 
when the stretching of the asperities is equal to x is given 
by N c P c (0 d£- The contacts that break when their 





FIG. 4: (Color online): The averaged final distribution Q s (x) 
for the parameters used in Figs.[2]and[3] The lower panel com- 
pares the functions P c (x) (solid line) and —dQ s (x)/dx (open 
diamonds). Red dotted curve and circles show the averaged 
final distribution of static thresholds. 



displacements are increased by AX are those which have 
their thresholds between x and x + AX, i.e. N C P C (X) AX. 
Therefore 



P(x) = 



Pc{x) 



(6) 



The broken contacts relax and have to be added to the 
distribution around x ~ 0, leading to the third term in 
Eq. (Q}. We denote by R(x) the normalized distribution 
of the stretching of the asperities that form a new con- 
tact after a slip event. Writing that all broken contacts 
described by AQ-(x; X) immediately reappear with the 
distribution R(x) (see Appendix |D|) . we get 



/oo 
dx' AQ_(x';X). 
-OO 



(7) 



Equation (0} can be rewritten as 



[Q(x; X + AX) — Q(x; X)] + [Q(x; X) - Q(x - AX; X)} 



= -AQ-(x;X) + AQ+(x;X). 

Taking the limit AX —> 0, we finally get the integro- 
differential equation 



dQ{x;X) dQ(x;X) 



dx 



dX 



P(x)Q(x;X) 



R(x) / dx' P{x')Q(x';X), 



(8) 
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which has to be solved with the initial condition 



C. Relation between the EQ and ME 



Q(x; 0) = Qini(a;). 



(9) 



Note that Qi n i(x) cannot be an arbitrary function, be- 
cause the contacts that exceed their stability threshold, 
must relax from the very beginning. 

In the following we select R(x) — 6(x), i.e. we as- 
sume that a broken contact sticks again only after a 
complete relaxation. This is physically reasonable and 
simplifies the analytical calculations, although this re- 
striction could be easily lifted for numerical solutions of 
Eq. ©. 

Integrating both sides of Eq. flSJ) over x from — oo to 
oo, we see that the normalization condition ([2]) is satis- 
fied, taking into account Q{— oo,X) — Q(+oo,X) = 
because contacts cannot be infinitely stretched without 
breaking. 

Once the distribution Q(x;X) is known, we can cal- 
culate the total force F with Eq. ((3]) and find the static 
friction force as the maximum of F(X), i.e., F s — F(X S ), 
where X s is a solution of the equation 



F'{X) = dF(X)/dX = 0. 



B. Steady state solution 



(10) 



The smooth-sliding solution Q s (x), i.e. the solution 
which does not depend on X, can easily be found directly 
[l4j . In the steady state Eq. ([5]) reduces to 



+ P(x)Q(x) = 5{x) dx' P{x')Q(x'). (11) 



dQ(x) 
dx 



In setting the lower bound of the integral we have as- 
sumed that P c (x) = for x < 0, which agrees with its 
physical meaning because, if x < 0, a positive variation 
AX actually reduces the absolute value of the force on a 
contact so that it does not cause its breaking. 

Its general solution can easily be derived for x > and 
leads to 



,(x) = C- 1 Q(x)E P (x), 



(12) 



where Q(x) is the Heaviside step function, O(x) =0 for 
x < and 6(x) = 1 for x > 0, 

E P {x) = e~ u <-*\ U(x) = [*dZP(Z) , (13) 

Jo 

and this solution also verifies the equation in the limit 
x — > 0. The normalization condition for Q s (x) gives 



C= I dxEp(x) 
lo 



The friction force is then equal to 

/>oo 

F = (N c (ki)/C) dxxEp(x). 



(14) 



(15) 



The EQ model is expressed in terms of the distribution 
of breaking thresholds of the contacts P c (x). Therefore, 
to connect the two approaches it is interesting to relate 
the steady state solution of the ME equation and the 
function P c (x). Using Eq. © we get 



= C -1 e(a:)exp 



Pc(0 



o .r luetic 



(16) 



The integration in the exponential can be read- 
ily performed since the integrand is of the form 

(d/df) In /"PcK'X- We § et 



i (x)=C- 1 Q(x) P c (0d£ 



or (x > 0) 



= - c -ip c (x) 

dx 



(17) 



(18) 



in agreement with the observations deduced from the nu- 
merical simulations of the EQ model (Fig.|4j lower panel). 
For the parameters used in Figs. ErEl the numerical con- 
stant is C = 0.9999. 

Notice also the useful relationships U'(x) = P(x), 
E' P (x) = -E P (x)P(x), E P (x) = 1 for x < 0, and 



P c {x) = P(x) E P (x) for x > 0, 



(19) 



Two simple examples are considered in Appendix [Bl 



the latter may be used in numerical solution of the ME 
instead of Eq. (j6]). 

For the simple example of a rectangular P c {x) distri- 
bution, which is considered in Appendix [B] the model 
admits an exact solution both for the EQ and ME ap- 
proaches. We checked that the earthquake model with 
the distribution (|B5[) and the master equation model with 
the expression (|B6[) for P(x) exactly have the same solu- 
tion for any initial configuration. 



D. Non-stationary solution 

The numerical solution of the master equation ([8]) 
when X starts to grow due to an external driving is pre- 
sented in Figs. [5] and [6] for a Gaussian distribution of 
static contact thresholds, P c (x) — G(x;x s ,(j s ). This 
solution may be compared with that of the EQ model 
of Figs. [2] and [3] for the same initial condition. One can 
see that they are almost identical, except for the noise 
on the earthquake model distributions. The distribution 
Q{x;X) always approaches the stationary distribution 
Q s (x) given by Eq. (|T2"|) . The final distributions of the 
EQ model and the ME approach are compared in Fig. [7] 

Let the distribution P c (x) be characterized by the av- 
erage value x s and the dispersion a s . Studying the evo- 
lution presented in Fig. [SJ we observe that every-time X 
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FIG. 5: (a) Solution of the master equation in the initial stage 
when X starts to grow. The model parameters are the same 
as in Fig. [2] (b) The corresponding dependence F(X). 
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FIG. 6: The same as in Fig. [5] for long times with the incre- 
ment AX = 1.09. 



increases by x s the distribution Q{x\ X) broadens so that 
its standard deviation grows by ~ a s . Therefore, any ini- 
tially peaked distribution tends to approach the station- 
ary one according to \Q(x;X) — Q a [x)\ cx exp(— X/X*), 
where X* ~ x 2 s /(X s . 

For one particular but important choice of the initial 
distribution, when all contacts are relaxed at the begin- 
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FIG. 7: (Color online): The final distribution Q(x) for the pa- 
rameters used in Fig. [B] (solid curve; crosses show the averaged 
final distribution for the earthquake model from Fig. [4| . The 
red dotted curve shows the distribution P c (x), blue broken 
curve - the function P(x). 



ning, Qi a i(x) = S(x), one can analytically express the 
initial evolution of the solution versus X. Namely, if 
P(x) — for the stretchings x < x*, then, as can be 
checked by its substitution in the master equation, the 
solution for the displacement X < 2x* is the following: at 
the beginning, for the displacement < X < x* , the solu- 
tion is trivial, Q(x;X) = 5{x-X) and F(X) = N c {h)X; 
for larger displacements, x* 



< X < 2x* c , the solution is 



Q(x;X) = 



where 



Q(x)b(X-x) 


a(X) S(X - x) 



for 0<x<X-x* cl 
for X — x* < x < x* c 
for x > x* 



b(0=P(0a(0 
and a'(£) = -P(0<»(0. or 



a(C) = exp 



C > <■ 



(20) 



(21) 



(22) 



Then one can calculate the friction force for the displace- 
ment interval x* < X < 2x* 



F{X) = N c {h) 



Xa(X) + (X-x* c ) 



x 



d£a(Z) 



(23) 

The solution (|2^|) allows us to find the static friction 
force F s — F(X S ) with Eq. (fit))) , which reduces to the 
equation X s a(X s )P(X s ) = 1 and corresponds to the first 
maximum on the variation of F(X). As was pointed out 
by Farkas et al. [HI], the ratio F s /Fk, where Fk is the 
kinetic friction, i.e. the plateau reached by F(X), takes 
values between 1 and 2 [larger values corresponding to 
narrower distributions P c (x)], and it is determined solely 
by the initial distribution Qi n i(x), reaching a maximum 
for Qi n i(x) = 5(x) (it was called "the initial coherence in 
strain distribution of the contacts" in Ref. Ilal). 
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Thus, in a general case, as X grows the solution of the 
master equation always approaches the smooth-sliding 
state given by Eqs. (fT2l - [T4|) . However, there is one ex- 
ception to this general scenario, when the model admits a 
periodic solution. This is the singular case when all con- 
tacts are identical, i.e., all contacts are characterized by 
the same static threshold x s , so that P c (x) = 8{x — x s ). 
The steady-state solution of Eq. © for this case is de- 
scribed in Appendix [E] We emphasize that the solution 
(|2T))) is valid for the continuous distribution P c {x) only, 
it cannot be used for the singular P c {x) — S(x — x s ) case. 



E. Nonrigid substrates: stick-slip 

The master equation allows us to compute the friction 
force F(X) when the bottom of the sliding block is dis- 
placed by X. But actually in an experiment one does 
not control X . The displacement is caused by a shearing 
force Ft applied to the top of the sliding block, which 
displaces its top surface by Xt- As the strain on the 
sliding block is usually small, its deformation can be as- 
sumed to be elastic, so that Xt is related to the applied 
force by 



F T =K(X t -X), 



(24) 



where K is the shear elastic constant of the solid block. 
The solution that we discussed above amounts to assum- 
ing that the sliding block was infinitely rigid, K — > oo, 
so that Xt = X at any time. In this case we found that 
the sliding always tends to a smooth-sliding steady state, 
as expected for a stiff system [H, H, 0| • Let us now con- 
sider the case of finite K < oo. The total force applied 
to the bottom of the sliding block, which determines its 
displacement X , is the sum of the applied force and the 
friction force 

F tot = K{X T -X)- F(X) . (25) 
It can be viewed as derived from the potential 
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x 



V(Xt, X) — —K(Xt — X) + F(®<%, (26) 

which determines the behavior of the sliding block sub- 
jected to friction and the applied force. 

A necessary condition for smooth sliding is that Xt 
and X grow together, with Xt — X = B, where B is 
a constant that measures the shear strain in the sliding 
block. This leads to the condition 



dV 



d(X 7 



X) 



dV 
dX 







(27) 



which simply means that the total force on the sliding 
interface vanishes. Smooth sliding also requires this state 
to be stable 



d 2 v 



d{X T - X) 



d 2 V 
dX 2 



> or F'{X) > -K . (28) 



These arguments exactly coincide with the "elastic in- 
stability" widely discussed in literature (e.g., see Ref. [3| 
and references therein). 

If we start from relaxed asperities, in the early stage 
of the motion F(X) is a growing function of X, and then 
it passes by a maximum when some contacts start to 
break and reform at lower asperity stress. Depending on 
the value of K, two situations are possible. For large K 
(stiff block), F' (X) never falls below ~K and the smooth 
sliding is a stable steady state. In this case the system 
evolves towards the regime F(X) = Const defined by 
Eq. ([15]). For small K (soft block), F'(X) can become 
smaller than —K and the stability condition ([28)) is no 
longer valid. In this case Eq. (|28|) . in the the limiting 
case of the equality, defines the maximal displacement 
X m that the contacts can sustain. A larger displace- 
ment breaks all the contacts simultaneously, causing a 
quick slip of the block before the contacts form again 
with a fresh distribution Q{x, X rn ) = R(x) and the same 
process can repeat again. The system may reach the 
regime of stick-slip periodic motion [l8| . Using the solu- 
tion (|23p for the delta-function initial distribution, we get 
F'(X) = N c (h) [1 - XP(X)a(X)}, so that we can find 
the period X ss of the stick-slip motion with Eq. (|28l) . 
Ignoring the slip time, it is given by the solution of the 
equation X a(X)P(X) = l + K/N c {h) for X > x* c . Note 
that in the case of an extremely soft substrate, K — > 0, 
the motion (almost always) corresponds to stick-slip, as 
it should for such a soft system [l|, y, 0] ■ 

Thus, depending on the distribution P c (x), the system 
demonstrates either stick-slip or smooth sliding. Smooth 
sliding is achieved if K > K*, otherwise the elastic 
instability occurs which may result in stick-slip. As 
K* = max[— F'(X)] and F(X) reaches its maximum for 
X ~ x s , one can estimate that K* ~ N c {k) x s /cr s ; thus, 
the ratio o s jx s controls the appearance of the elastic in- 
stability. 



IV. TEMPERATURE EFFECTS 

The effect of a nonzero temperature is connected with a 
change of the fraction density of breaking contacts P(x) 
in the master equation flSJ as first discussed by Pers- 
son [l4j. Indeed, for a single contact with the static 
threshold a: s at zero temperature, the contact does not 
break at all for x < x s . But when T > 0, the contact 
may relax due to a thermally activated jump before the 
threshold x s is reached. The rate h(x;x s ) of this process 
is defined by 



dQ(t)/dt = h(x;x s )Q(t), x< 



(29) 



where Q(t) is the probability that a contact existing at 
t = is not thermally broken at time t. For a set of 
contacts, Eq. (f2T)]) has to be generalized to 



dQ(t)/dt = H(x)Q(t) 



(30) 
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with 



H(x) 



dx h(x; x ) P c (x'). 



(31) 



For a sliding at velocity v so that X — vt, the thermally 
activated jumps can be incorporated in the master equa- 
tion, if we use a corrected expression Pt{x) defined by 



P T {x) = P(x)+H(x)/v 



(32) 



instead of the zero-temperature breaking fraction density 

The rate of thermal activation of the contacts, h{x; x s ) 
in Eq. (f2l?|). can be estimated by the Kramers relation. 
Let V(x) be the binding energy of the contact, and 
AE(x; x s ) — V(x s ) — V(x), the activation energy for the 
contact to break. For "soft" , or "weak" contacts, when 
AE(0;x s ) > k B T, the rate h(x;x s ) is given bv [T3. fl5l| 



h{x\ x s ) ~ u exp [— AE(x; x s )/kBT] , 



(33) 



where u is a prefactor corresponding to the attempt fre- 
quency. For an overdamped dynamics of the contacts 
u « ujq/2tti] with the characteristic frequency uj 2 ~ 
k/m ~ c 2 1 A (m is the contact mass) and the damping co- 
efficient rj ~ c/y/A which gives ui ~ c/2tt^/A ~ 10 10 s _1 
[Tij . If we assume that the binding potential of a con- 
tact can be approximated by the elastic properties of this 
contact, then V(x) = \ kx 2 . In this case the activation 
energy for contact breaking takes the form [lij 



1 



AE(x;x s ) = ^k(x 2 



x ) ~ kx s (x s x), 
and the function H(x) in Eq. (|32[) is given by 



H(x) 



to e 



kx 2 /2k B T 



-k(, 2 /2k B T 



(34) 



(35) 



For "stiff", or "strong" contacts, when AE(0;x s ) > 
fcflT, contact breaking only occurs with a significant 
probability when the stretching x of a contact is close to 
the threshold x s . The harmonic expression of the binding 
energy can no longer be used and a cubic approximation 
is more appropriate. This leads to a correction in the 
activation energy [22|, [23[ 

AE(x; x s ) = A£(0; x B ) (1 - x/x s f 2 , (36) 

and a renormalization of the prefactor uj, 

u lj (1 - x/x s ) 1/2 , (37) 

so that the function H(x) is given by 

, , f°° x\^ kfa-x/o*' 

(38)" 

The contributions (|33|) or (|3"5)l to the rate Pt{x) lead 
to appearance of a low- a; tail. Its magnitude grows with 
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FIG. 8: (Color online): The rate Pt{x) for soft contacts at 
different temperatures T = (black short dashed), 0.03 (ma- 
genta dash-dotted), 0.1 (blue dotted), 0.3 (red dashed) and 1 
(black solid curve) for a given velocity, u/kv — 1. Inset: 
the rate Pt(x) at different velocities uj/kv — (black short 
dashed), 0.3 (magenta dash-dotted), 1 (blue dotted), 3 (red 
dashed) and 10 (black solid curve) for a fixed temperature 
T = 0.3. The red short-dotted curve shows the distribution 
P c (x) (Gaussian with x s — 1 and a 3 = 0.05). 
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FIG. 9: (Color online): The steady-state distribution Q s (x) 
for T = 0.3 and different velocities uj/kv = 0.3 (magenta dash- 
dotted), 1 (blue dashed), 3 (red short-dotted) and 10 (black 
solid curve) for the soft contacts. The distribution P c (x) is 
Gaussian with x s = 1 and a 3 = 0.05. 



temperature as well as with a decrease of the driving 
velocity v as demonstrated in Fig. [8] 

Raising temperature leads to a shift of the distribution 
Q(x) to lower values, so that the friction force decreases 
when T grows. This effect is larger for lower sliding veloc- 
ities as shown in Fig. [HI In the limit v — > 0, all contacts 
finally break if T ^ 0, so that Q s (x) — !• S(x) and F fc -» 0; 
in this limit we have a smooth sliding associated to a 
creep motion of the contacts. 

The variation of Fk(T) versus T for different sliding 
velocities in the smooth sliding regime for soft and stiff 
contacts is presented in Fig. [TUh . As expected, the fric- 
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FIG. 10: (Color online): (a) Dependence of the kinetic friction 
force Fk in the steady state on the temperature T for different 
velocities v = 0.01 (black squares), 0.1 (red circles), 1 (ma- 
genta up triangles) and 10 (blue down triangles), (b) Fk(v) 
for different temperatures T = 0.01 (black squares), 0.1 (red 
circles), 0.3 (magenta up triangles) and 1 (blue down trian- 
gles); the inset shows the same in log- log scale. Solid curves 
and symbols are for soft contacts, and broken curves and open 
symbols, for stiff contacts. The distribution P c (x) is Gaussian 
with x a = 1 and a 3 = 0.05; ui/k — 1. 



tion force decreases when T grows and tends to zero when 
T — » oo. Figure [TOb shows the dependence of kinetic fric- 
tion on the sliding velocity. The force Fk monotonically 
increases with v, approaching the T = limit when v — ¥ 
oc. The dependence Fk(v) at smooth sliding was esti- 
mated analytically by Persson [3] (see also Ref. (24[) who 
found Fk(v) oc v in the low- velocity limit, Fk(v) oc ]nv for 
intermediate velocities, and Fk(v) — F(oo) oc —ujT 2 jkv 
in the high-velocity case. 

Thus, for T > the friction force Fk increases 
with the driving velocity, i.e. dFk(v)/dv > 0, which 
stabilizes the smooth sliding regime. However, ex- 
periments show that the temperature-induced velocity- 
strengthening only dominates the aging-induced velocity- 
softening (see Sec. E} at high velocities v > 10 — 
100 /xm/s Q. 

Of course, a nonzero temperature also affects the dy- 
namics of the approach to the steady state as shown in 
Fig. [TT] for different driving velocities. The higher the 
temperature and/or the lower the velocity are, the lower 
is the static friction force F s [determined by the first max- 
imum of the F(X) dependence], and the faster F(X) 
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FIG. 11: The dependences F(X) for T = 0.3 at different driv- 
ing velocities v, uj/kv — 0, 0.3, 1, 3 and 10. The distribution 
Pc(x) is Gaussian with x s = 1 and a s = 0.05, the initial dis- 
tribution Qini(a;) is Gaussian with xi n i = 0.1 and <Ti n i = 0.025. 



approaches the steady-state smooth sliding. It is impor- 
tant to consider the first cycle of the F(X) dependence, 
which defines the lowest value of F'(X) because it de- 
termines the stability limit of the smooth sliding regime 
according to Eq. fl28} . As shown in Fig. higher tem- 
peratures (Fig. [T2h ) raise the minimum value of F'(X), 
which in turn extends the interval of the model param- 
eters for which the smooth sliding takes place. At the 
same time, the higher is driving velocity (Fig. H2b). the 
smaller is the minimum of F'(X), so that the narrower 
is the interval of model parameters where the smooth 
sliding regime exists. Thus, we come to a surprising con- 
clusion that, at T > 0, reducing the pulling velocity may 
lead to a transition from stick-slip to smooth sliding or 
creep motion, while one generally expects that smooth 
sliding is reached when the velocity increases. However 
a transition to smooth motion by reducing velocity has 
also been observed experimentally 25j . 

The behavior of Fk on T and v described above is only 
qualitative. A quantitative study to identify the temper- 
ature range in which these effects would play a signifi- 
cant role in an actual material would require a specific 
study. However it qualitatively agrees with tip-based ex- 
periments [Ifl which suggest that, in these exper- 
iments, the tip/substrate contact occurs through many 
atoms, i.e. it actually corresponds to multiple contacts. 
Surprisingly, the dependences obtained within the ME 
approach, perfectly agree with the experimental ones for 
a model tribological system [281 - 131) . where the wearless 
kinetic friction of the driven lattice of quantized mag- 
netic vortexes in high-temperature cuprate superconduc- 
tors was studied (e.g., compare inset of Fig. [TUb with 
Fig. 3 of Ref. [H|). 

However our approach only claims to explain the gen- 
eral trends of the behavior of tribological systems, be- 
cause we neglected several important factors. First, we 
neglected contact aging (see Sec. |V| . Second, we as- 
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FIG. 12: (Color online): The same as in Fig. El for short X. 

(a) is for different temperatures T — (red dash-dot-dotted), 
0.03 (magenta dash-dotted), 0.1 (blue dashed), 0.3 (red dot- 
ted) and 1 (black solid curve) at the fixed velocity, oj/kv = 1; 

(b) is for different velocities v, ui/kv — (red dash-dot- 
dotted), 0.3 (magenta dash-dotted), 1 (blue dashed), 3 (red 
dotted) and 10 (black solid curve) at the fixed temperature 
T = 0.3. 



sumed that R(x) = S(x), while at T > the distribu- 
tion R(x) of stretchings when contact form again after a 
slip event should correspond to the Boltzmann distribu- 
tion around x — with a width depending on T [ID, [24| • 
Third, in a real system the contacts are heated due to 
sliding and even may change their structure. Forth, we 
ignored the elastic interaction between the contacts. Fi- 
nally, in macroscopic experiments the wearing of sub- 
strates in sliding contact may mask the temperature and 
velocity dependences. A more detailed investigation of 
these aspects deserves separate studies. 



AGING OF THE CONTACTS 



In the approach described in Sec. IIII1 the function 
F(X) was independent on the driving velocity because 
we neglected two time-dependent effects: first a broken 
contact needs a minimum delay r m to stick again, sec- 
ond we ignored the aging of the contacts once they are 
formed, which contradicts most of the experimental re- 
sults [J H, El ■ Experiments as well as molecular dynamics 
simulations indicate that the static friction force grows 
with the lifetime t w of stationary contacts, i.e. the time 
interval during which they stay static prior to sliding. 

The simplest way to obtain a velocity dependence of 



the function F(X) is to take into account the delay time 
T m . This is considered in Appendix [Dj In the most 
general case the delay time may not be the same for all 
contacts but exhibits some distribution. In Ref. [IH we 
showed that the condition r m > is the second necessary 
condition for the appearance of stick-slip. 

Let us now consider the aging of contacts which re- 
quires a more involved analysis because we must take into 
account a time evolution of the distribution P c (x) of the 
contact breaking thresholds. Let the newborn contacts 
be characterized by a distribution P C i(x) with the aver- 
age value x a i and the dispersion cr s ;. If aging effects are 
taken into account, then typically x s grows with time, 
while a s decreases, because the area of contact slowly 
grows under pressure. As it does it faster for the small 
contacts, which are less resilient, this tends to reduce the 
dispersion in the properties of the contacts. As a result, 
at t — > oo the distribution P c (x) approaches a distribu- 
tion P c f(x) with x s f > x s i and a s f < a S i. The final 
distribution is a property of the material. It is there- 
fore a known quantity for the model once the material 
has been specified. If we assume that the evolution of 
P c (x) corresponds to a stochastic process, then it should 
be described by a Smoluchowsky equation 

dP c - d ( , x d 

-^=DL x P Cl where L x = — lB(x) + — 

" (39) 

in which the "diffusion" parameter D describes the rate 
of aging, B(x) = dU(x)/dx, and the "potential" U(x) de- 
fines the final distribution, P c f(x) oc exp —U{x) ; there- 
fore we can write 



B(x) 



dP c f (x) I dx 

Pcfix) 



(40) 



Then, the equation dP c /dt = D L X P C naturally leads 
to a growth of the average static threshold from x S i to 
x s f with the time of stationary contact, as widely as- 
sumed in earthquake-like models of friction 0, @, EH- 
A characteristic timescale of aging may be estimated as 

Paging — K^sf •Esi) I D . 

At the same time, the contacts continuously break 
and form again when the substrate moves, as described 
by third and forth terms in Eq. ([§]). This introduces 
two extra contributions in the equation determining 
8P c (x;X)/dX in addition to the pure aging effect de- 
scribed by Eq. (|3"9")l : a term P(x; X) Q(x; X) takes into 
account the contacts that break, while their reappearance 
with the threshold distribution P C i{x) gives rise to the 
second extra term in the equation. Combining the equa- 
tion describing the evolution of the distribution Q(x;X) 
with the equation describing the aging of the contact 
properties, and taking into account the relation (TT91 . we 
come to the system of three equations: 



dQ{x-X) dQ{x;X) 



dx 



dX 



P{x;X)Q(x;X) 
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FIG. 13: Evolution of P c {x) due to two concurrent pro- 
cesses: the first process is the aging of asperities from the 
initial (fresh) distribution P C i(x) (Gaussian with x s i = 0.5 
and a si = 0.05) to the final distribution P c f(x) (Gaussian 
with x s f = 1 and a s f = 0.02), and the second one is the 
breaking/re-binding of the contacts. Dx = D/v = 5 x 10 -4 , 
the initial distribution Q- m i(x) is Gaussian with Xini = 0.1 
and o"ini = 0.025. The curves are plotted with the increment 
AX = 0.05. 



= 6(x) I dx' P(x';X)Q(x';X), 

dP c {x-x) 



dX 



D x L x P c {x; X) + P(x; X) Q{x; X) 



P cl (x) / dx' P(x';X)Q(x';X), 



P c (x;X) = P(x;X) exp 



Jo 



(41) 



where D x = D/v, and v = dX[t)/dt is the driving veloc- 
ity. Thus, in the case of high driving (v — > oo, Dx — > 0) 
we recover the previous behavior with P c {x) = P c i(x). In 
the case of low driving (v — > 0, Dx —> oo) we again ob- 
serve the same type of behavior but with P c (x) = P c f(x), 
and in the case of < Dx < oo we have an interplay 
of two processes, the aging which tends to change P c {x) 
towards P c f (x) , and the breaking- reformation of the con- 
tacts which returns P c (x) to P C i(x). This is illustrated 
in Fig. [T31 where we plot the evolution of P c (x) when X 
continuously increases. As a result, the friction force F 
has to depend on the sliding velocity as shown in inset 
of Fig. IT4"b. Because typically x S i < x s f, the kinetic fric- 
tion force Fk generally decreases when v grows, so that 
dF k (v)/dv < 0, which may lead to an instability of the 
smooth sliding regime. 

The steady-state solution of Eqs. (l4"Tj) can be found 
analytically in the low- and high-velocity cases. In the 
v — > limit, 



P c (x) « P cf (x) 
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FIG. 14: (Color online): (a) The final distribution Q s (x) for 
five values of the parameter Dx = D/v: 10~ (black circles), 
10~ 4 (blue up triangles), 3 x 10 -4 (red crosses), 10~ 3 (ma- 
genta down triangles), and 10~ 2 (black diamonds), (b) The 
corresponding distributions P c (x). Larger open circles show 
the initial distribution P C i(x) (red dashed) and the final dis- 
tribution Pcf(x) (blue dotted). The inset displays the depen- 
dence of the kinetic friction force Fk in the steady state on 
the driving velocity v = 1)- Pd{x) is Gaussian with 

x si = 0.5 and a s i = 0.1, and P c f(x) is Gaussian with x s f = 1 
and a s f = 0.01. 



jf <% p-/(o £ [p cf (e) ~ p ci (a , 

which leads to Fk(v) — F k (0) oc —v/D, while in the high- 
velocity case, v — > oo, we have 



P c (x) ^ P cl {x) + - C[Pi}^- 
v dx 



B(x) 



d 
dx 



Pci(x), 



which gives F k (v) — Fk(oo) oc D/v. In these expres- 
sions C[P] designates the constant C defined by Eq. (fl4| . 
which takes different values depending on the expression 
of P(x). Both limiting cases may be combined within 
one fitting formula, 



F k (v) «F fe (oo) + 



F fc (0)-F fc (oo) 

l + t'/Waging 



(42) 



where Waging defines the velocity below which aging effects 
become essential. For example, for the parameters used 
in Fig. 1141 we found that f aging ~ 6.5 x 10 3 L>. 

When the sliding substrates are not rigid, the effect of 
contacts aging may lead to a transition from stick-slip to 
smooth sliding with the increase of the sliding velocity. 
When the driving velocity increases and reaches a criti- 
cal value v = v c the transition is abrupt (first-order) as 
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FIG. 15: The dependences F{X) for two driving velocities 
v — 1.5 (a) and v = 1.6 (b). The initial (fresh) distribution of 
thresholds P C i(x) is Gaussian with x S i = 0.5 and a a i = 0.05, 
the final distribution P c f(x) is Gaussian with x s f = 1 and 
cr s f = 0.02. The initial distribution for contacts Qi n i(x) is 
Gaussian with X[ n i = 0.02 and <7i n i = 0.01 (for calculation 
reasons, the same is used for the "refreshed" distribution after 
the break at the point F'(X) = —K). The parameters are 
(ki) = 1, K = 3.5N c (ki) , and the rate of aging is determined 
by the coefficient Dx = D/v with fl = 5x 10~ 4 . 



contact thresholds, which may be time dependent as dis- 
cussed in Sec. [V] on aging. Let us now discuss possible 
physical origins of this distribution. 

To establish the master equation, it is convenient to 
express the properties of the contacts in term of their 
stretching x, hence P c {x), but at the microscopic scale 
contacts are generally characterized by their breaking 
force. The distribution P c (x) can be related to the dis- 
tribution P c (fs) of the static friction force thresholds of 
the contacts. If a given contact has an area A, then it is 
characterized by the static friction threshold f s oc A and 
the elastic constant k ~ p <?\f~A, where p is the mass den- 
sity and c the sound velocity (assuming that the linear 
size of the contact and its height are of the same order 
of magnitude [3). The displacement threshold for the 
given contact is x s = f s /k, so that 

/,«a& (43) 

or df s /dx s oc x s . Then, using the relationship 
P c {x s ) dx s = P c (fs) df s , we obtain 

P c (x s ) oc x s P c [f s (x s )]. (44) 
A. Rough surfaces: the GW model 



demonstrated in Fig. [15j Moreover, the system exhibits 
hysteresis: starting from the smooth-sliding regime, if the 
velocity decreases we expect a smooth sliding to survive 
down to velocities much lower than v c . 

In a general case, the parameter D in Eq. (jTTj) 
and, therefore the rate of contacts aging, depends on 
the system temperature T. Typically, the aging rate 
should increase with T according to Arrhenius law, 
D oc exp (— e/ksT), where e is an activation energy for 
this process. Besides, aging mechanisms may be much 
more involved than the simplest one described by the 
Smoluchowsky equation (|39[) and may even consist of 
two stages [1] , the "geometrical" growth of contacts (de- 
scribed, e.g., by the Lifshitz-Slozov theory, see Sec. lVIE|) 
and then a structural reordering of asperities. 



VI. ORIGIN OF THRESHOLD DISTRIBUTION 

As discussed in the introduction the formalism using 
the master equation to describe friction allows us to sep- 
arate the calculation of the friction force, assuming that 
the statistical properties of the contacts are known, from 
the analysis of the properties of the contact themselves. 
In the previous sections we focused on the calculation of 
the friction force, which is our main goal. However it 
is interesting to examine the properties of the contacts 
because it allows us to draw some general conclusions on 
friction, particularly on the possible existence of stick- 
slip. Within our approach the properties of the contacts 
are entirely described within the distribution P c (x) of the 



Let us consider the contact of two nominally flat sur- 
faces, for which the appararent contact area is large so 
that individual contacts are dispersed and the forces act- 
ing through neighboring spots do not influence each other 
21|. Their study can be cast into the problem of the 
contact of a rigid plane and a composite surface whose 
topography is the sum of both topographies, with an 
appropriate renormalization of some parameters such as 
the Young modulus, so that the contact problem reduces 
to the statistics of independent asperities [4(. Let the 
rough surface be characterized by hills of heights {hi} dis- 
tributed with a probability Ph(h). Following the Green- 
wood and Williamson (GW) model of the interface [2lj |. 
let us assume that all hills have a spherical shape of the 
same radius of curvature r. When this surface is pressed 
with another rigid flat surface, which takes a position at 
the level ho, then the hills of heights h > ho will form 
contacts, or asperities. If the contacts are elastic (the so 
called Hertz contacts [32]), then the contact of height h 
has the compression (h — ho), its area is irr(h-ho), and it 
bears the normal force fi(h) « (4 7 r/3)£;*r 1 / 2 (/i - h f' 2 , 
where E* is the effective Young modulus (if both con- 
tacting substrates have the same Young modulus E, then 
E* = £7/2(1 - v 2 ), where v is the Poisson modulus [32|). 
It is reasonable to assume that the shear static threshold 
for the given contact is proportional to the load force, 
f s (h) « pfi(h), or 

/ s = (4 7 r/3) M £;*r 1 /2(/ l _/ l0 )3/2 ) (45) 

where p < 1 is a constant. Then the distribution of 
static thresholds P c (f s ) can be coupled with the distribu- 
tion of asperities heights Ph{h) with relation P c (f s ) dfs oc 
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P h (h) dh, or P (J„) cx {dh/df s )P h [h(f s )], where dh/df s <x 
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FIG. 16: (Color online): Dependence of the friction force F 
on X for elastic (£?' = 1.588, solid curve) and plastic (B" = 
1.251, broken curve) contacts with exponential distribution 
of heights. The initial distribution Qi n i(x) is Gaussian with 
xi n i = and aim = 0.025. The inset shows the corresponding 
distributions P c {x), Eqs. ([47j) and (|48)l . 
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FIG. 17: (Color online): Dependence of the kinetic friction Fk 
in the smooth sliding state on the parameter B for the elastic 
(solid curve) and plastic (broken curve) contacts. Inset: the 
same in log-log scale; dotted lines show power-law fits. 

For a strong load, when the local stress exceeds the 
yield threshold Y, the contacts begin to deform plas- 
tically. When all contacts are plastic, the local pres- 
sure on contacts is pi oa d = H, where H is the hardness 
[H k, 3Y for the spherical geometry of asperities; for 
metals H « (1CT 3 - IQ- 2 )E - 10 9 Pa]. Then the normal 
force at the contact is fi(h) w nr(h — ho)H. Assuming 
again that f s (h) « jifi (h) , we obtain 



f s = nr(h - h )nH, 
so that P c (fs) oc P h [h(f s )] now. 



(46) 



For example, if we consider the exponential dis- 
tribution of heights introduced by Greenwood and 
Williamson J2l| to get analytical results Ph(h) — 
hr 1 exp(—h/h) 0(h), where h is the average height, then 
with Eq. (|44|) for the elastic contacts we obtain (/, x > 0) 



Pc(f)*r 1/3 eM-Bif 2/3 ), or 
P c {x) cx a; 1/3 exp(-B'a; 4/3 ), 



(47) 



where B l = [(4 7 r/3) 2 / 3 /ir 1 /3(^*)2/3] 
plastic contacts 



while for the 



P c (f) oc &q>(-B 2 f), or P c (x) cx x exp(-B"x 2 ), 

(48) 

where B>2 = (irhrfiH) . B' and B" are numerical con- 
stants depending on the material. 

The distributions P c (x) for elastic and plastic contacts, 
Eqs. (|4"T)) and (|4"5)) . are presented in Fig. [TTdI (inset). We 
see that the force F(X) increases (almost) monotonically 
with X, approaching the kinetic value F/.. Thus, in the 
case of a contact of rough surfaces, a relatively high con- 
centration of low-threshold contacts prevents the stick- 
slip motion even for a very soft tribosystem (i.e., even 
in the case of a very low elastic constant of the driv- 
ing spring K). The kinetic friction force F/- depends 
on the parameter B according to the power law (see 
Fig. [TTJ), F fc oc B~y 4 oc h 3/4 for the elastic contacts, 
and Fk oc B^ 1 !" 1 oc h 1 / 2 for the plastic contacts. 

Similarly one can find P c (f) for a more realistic Gaus- 
sian distribution of heights of the asperities [2l[ . In this 
case, the peak in the P c (x) dependence becomes nar- 
rower and higher, and the stick-slip regime may exist. 
Note that typical values for rough surfaces are of the or- 
der r ~ 10 — 100 /im, h ~ 1 /im, so that the average size 
of the contact is a w (rft,) 1 ^ 2 ~ 3 — 10 

Whether the contacts are in the plastic or elastic 
regime, depends on the dimensionless parameter ip = 
(E/Y)(h/r) 1 ^ 2 : ^> 1, which is typical for metals, cor- 
responds to the plastic regime while ip < 1 leads to an 
elastic regime, as, for instance, for rubber friction. Poly- 
meric glasses belong to an intermediate situation, -0^1, 
where only a fraction of contacts is plastic H. 



B. Rough surfaces: Persson's model 

The GW model used above completely ignores the elas- 
tic interaction between the contacts thus overestimating 
the role of low-threshold contacts. Recently Persson [33[ 
developed a contact mechanics theory which (indirectly) 
includes elastic interactions and leads to the correct low- 
threshold limit P c (f — > 0) = 0. Although Persson's the- 
ory is only rigorous for the case of complete contact (e.g., 
contact of rubber surfaces at high load) , it leads to a very 
good agreement with experiments and simulations even 
for a contact of stiff surfaces and low loads [34| • Not going 
into details, note that the distribution of normal stresses 
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er„ at the interface can approximately be described by an 
expression (a n > 0) 



in(fn) oc exp 



P~n — Q~n 
Acr„ 



-exp 



Acr„ 



(49) 

where <7„ is the nominal squeezing pressure, the distribu- 
tion width is given by Aa n = E*1Z X / 2 , and the parameter 
1Z is determined by the roughness of the contacting sur- 
faces, 



U= -L / dqq 3 I d 2 x (//(x)/((0)}( - ,c ' x 



(50) 



Assuming that a local shear threshold is directly propor- 
tional to the local normal stress, f s (xcr n , and again using 
Eq. (|4"4"|) . we finally obtain 



P c (x) oc xi exp 



2 -2 \ 2 

ar — ar 



exp 



2 , —2 \ 2 
2T + 0T 



(51) 



where a oc (E 1 *) 1 / 2 ^ 1 / 4 . The distribution flSJ) is char- 
acterized by a low concentration of small shear thresh- 
olds, P c (x) oc x 3 at x — >• 0, and a fast decaying tail, 
P c (x) oc exp[— (cc/cr) 4 ] at x — > oo, i.e., the peaked struc- 
ture of the distribution (|5ip is much more pronounced 
than in the GW model, Eqs. (gTJ) and (l48l) . 

The numerical results are presented in Fig. [18] for two 
values of the roughness parameter. In contrast to the 
results of the GW model, they show a non-monotonous 
behavior of F(X) which is compatible with a stick-slip 
motion. 




FIG. 18: (Color online): Dependence of the friction force F 
on X for the Persson model with x = 1, a = 0.5 (solid) 
and a = 1 (broken curve). The initial distribution Qi ni(^r) is 
Gaussian with X[ n i = and a ln i — 0.025. The inset shows the 
corresponding distributions P c (x), Eq. (|5ip. 



C. Flat surfaces: dry friction 

Let us now consider the dry contact of two flat surfaces. 
If both surfaces have an ideal crystalline structure, we 
get the singular case discussed in Appendix [E] However 
such a situation is exceptional. A real surface always 
consists of domains, which are characterized by different 
crystalline orientation or even different structures. 

MD simulations show a large variation of the static 
friction with the relative orientation of the two bare sub- 
strates [H, HH , although surface irregularities as well as 
fluctuations of atomic positions at nonzero temperatures 
may mask this dependence !37]. A variation of the fric- 
tion with the misfit angle was also observed experimen- 
tally, for instance in the FFM experiment made by Di- 
enwiebel et al. 38]. 

In order to estimate a possible shape of the function 
P c {x) resulting from the domain structure of substrates, 
let us consider a rigid island (domain) of size N a with 
a triangular lattice placed onto a 2D periodic potential 
with square symmetry V(x,y) = since + siny created by 
the bottom substrate. The atomic coordinates of domain 
atoms are Xi = X + (i + j/2)a and y j = Y + jh, where 
h = ay/3/2, i = -j/2, -j/2+m-l, j = 0, . . .,n ~l, 
riirij = N a (nja ~ rijh), and X, Y are the center of mass 
coordinates. If the island is rotated by an angle <fi, then 
the atomic coordinates change to x\ = Xi cos 4> — jji sin <p 
and yi = Xi sin <fi + jji cos 4>- For a fixed misfit angle 0, 
the total potential energy of the domain is U(X, Y) = 
^2ijV{xi,yj). Extrema (minima and saddle points) of 
the function U(X, Y) arc defined by dU/dX = dU/dY = 
0; then the activation energy for domain motion is given 
by e a = Laddie - U min . Assuming that f s - s a /a s oc e a , 
we can estimate the threshold force f s as a function of the 
misfit angle <f> and the domain size N a . This protocol was 
realized by Manini and Braun [3!| . The calculation of a 
histogram of the function e a (4>) leads to the distribution 
Pc(fs), assuming that all domains have the same size 
N a and that all angles are equally present. Next, one 
may average over different domain sizes N a , e.g., with 
a weight function w(N a ) — e~ Na ' Na , where N a is the 
average domain size. The distribution P c (f) obtained in 
this way may be crudely approximated by the function 
P c (f) oc exp(- ///). Assuming that Eq. (g?]) holds, we 
come to the distribution 



P c (x) oc x exp(— x 2 /a 2 ) , 



(52) 



which is similar to that of Fig. [TC] characteristic for 
rough surfaces. Therefore in this system we again expect 
smooth sliding without stick-slip due to a large concen- 
tration of configurations with small barriers. 

However, in the estimation we supposed that all an- 
gles are present with equal probability, while some an- 
gles could have preference due to their lower potential en- 
ergy. Calculation shows that this would suppress the low- 
barrier thresholds so that the distribution P c {x) would 
become more peaked, and the stick-slip could exist. 



15 



D. Lubricated surfaces: MD simulation 



E. Lubricated surfaces: Lifshitz-Slozov coalescence 



The dry-friction considered above is exceptional. In a 
real system, there is almost always a lubricant between 
the sliding surfaces (called "the third bodies" by tribol- 
ogists) which is either a specially chosen lubricant film, 
grease (oil), dust, wear debris produced by sliding, or wa- 
ter or/and a thin layer of hydrocarbons adsorbed from 
air. 

There is a large number of experimental and simula- 
tion studies of lubricated friction (e.g., see Refs. [HQ and 
references therein) . A thin confined film (less than about 
six molecular diameters of thickness) typically solidifies 
thus producing nonzero static friction. The value of f s 
strongly depends on the film thickness and its structure, 
and moreover it may change with the time of station- 
ary contact. In particular, Jabbarzadeh et al. [40l | have 
done the MD simulation of a thin lubrication film of do- 
decane. The transition from bulk liquid to high viscosity 
state when thickness decreases occurs at the thickness 
of six lubricant layers, and it appears due to formation 
of isolated crystalline bridges between the mica surfaces 
(across the film). As the thickness decreases further, 
these bridges increase in number and organize themselves 
into a mosaic structure with a long range orientational 
order. However in the previous MD study of the six-layer 
dodecane film between the mica substrates Jabbarzadeh 
et al. 5l| found also a "layer-over-layer" sliding regime 
with a very low friction. Such a configuration is found to 
be thermodynamically stable contrary to the metastable 
"bridge" configuration mentioned above. Thus, if the 
film is not ideally homogeneous along the interface but 
consists of domains of different structures and may be 
different thicknesses, it should be characterized by a dis- 
tribution of static thresholds. 

Moreover the top and bottom surfaces may be misori- 
ented as well, as was mentioned above. In particular, He 
and Robbins studied the dependence of the static [42[ 
and kinetic [43| friction on the rotation angle of the sub- 
strates for a lubricated system. It was observed that 
the static friction exhibits a peak at the commensurate 
angle (0 = 0) and then is approximately constant, the 
peak/plateau ratio being about 7 (for the monolayer lu- 
bricant film, for which the variation is the strongest). 

A detailed study of the dependence of the static fric- 
tion on the rotation angle <fi for lubricant films of thick- 
ness from one to five layers has been done by Braun and 
Manini [44j . The value of / s varies with <fi by more than 
one order in magnitude. When averaged over (f> (assum- 
ing that all angles are present with the same probability) 
and also over film thickness with some weighting factor, 
the distribution of static thresholds may approximately 
be described by the function P c (f) oc /exp(— ///). As- 
suming again that Eq. (|44p holds, we come to the distri- 
bution 



P c (x) oc x 3 exp(-x 2 /u 2 ) 
which would allow stick-slip. 



(53) 
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FIG. 19: (Color online): Evolution of the distribution P c (x) 
with the time of stationary contact for the Lifshitz-Slozov 
coalescence mechanism (a = 1 and B — 1). 

In the conventional melting-freezing mechanism of fric- 
tion, the lubricant is melted during slip, and solidi- 
fies when the motion stops. The solidification process 
can be described by the Lifshitz-Slozov theory (see [ijj . 
Sec. 100). At the beginning, grains of solid phase emerge 
within the liquid lubricant film. Then the grains grow in 
size according to the law f oc i 1 / 3 , where f is the average 
grain radius. The distribution of grains sizes is described 
by the following function: the number of grains with the 
radius from r to r + A?' is equal to PlsI 7 "/^) Ar/f, where 
the function Pls(u) is zero for u > 3/2 (so that the max- 
imum size of the grains is 1.5f), while lower sizes are 
distributed as 



Pls(u) oc 



u 2 exp [ 



-1/(1 -2u/3)] 



(u + 3) 7 /3(|_ M ) 



11/3 



(54) 



Due to coalescence of grains, the total number of grains 
decreases with time as N(t) oc i -1 . When the size 
of a grain exceeds the film thickness d (that will oc- 
cur when f{t) > d/3 after the delay time r m oc d 3 ), 
such a grain will pin the surfaces. Using the relation- 
ship P c (f)df oc PLs(r/f)dr/f, we obtain that P c (f) oc 
(dr /df)Phs{ r /r) / f. Because one grain gives the static 
threshold proportional to the contacting area, f s oc 
7r(r 2 - d 2 /4) for r > d/2, we have dr/df oc / _1 / 2 . Com- 
bining all things together, we obtain that 



P c (x) = P LS (u), u = p- 1 (l + Bx 2 



,1/2 



(55) 



where p{t) = 2f(t)/d — at 1 ' 3 (the pinning begins when 
p > 2/3), and B is determined by the system parameters 
(elasticity of the contacts, proportionality between the 
threshold f s and the area of the contacting grain, and the 
thickness of the lubricant film) . The distribution ([5"5j) is 
shown in Fig. Q1JJ Its shape suggests that it should lead 
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to the conventional tribological behavior: the stick-slip 
motion at low driving velocity and smooth sliding at high 
velocities. 

Thus, the described mechanism leads to a non-singular 
distribution of static thresholds and the typical stick-slip 
to smooth sliding behavior even for the case of ideal flat 
surfaces such as, e.g., the mica surfaces used in SFA ex- 
periments. Note that it incorporates both the delay time 
discussed in Appendix [D] and the aging of contacts con- 
sidered in Sec. [Vj the latter, however, is a more com- 
plicated process than follows from the simplest Smolu- 
chowsky theory. 

In a general case, we also have to take into account 
that the lubricant film may consist of grains (domains) 
of different orientations or even different structures. In- 
deed, the proportionality coefficient in the relation f s oc 
7r(r 2 — d 2 /4) used above, should depend on the misfit 
angle between the lubricant domain and the substrate, 
so that the distribution Pls( u ) introduced above, addi- 
tionally should depend on the misfit angle </>, Pls(m;</>). 
Thus, if there are grains with different orientations, dis- 
tributed according to a function i?^(0), then 

P c (x) = f d<t>R4,(<l))P hS (u;4>). (56) 

For example, if there are only two energetically equal 
possible orientations with (j) — and cf> — 7r/4, then 
i**(0) = 5«(0) + 5*fa-f/4). 

VII. CONCLUSION 

We introduced and analyzed in detail the earthquake- 
like model with a distribution of static thresholds. Re- 
ducing the problem to a master equation, we were able 
to find the exact solution numerically and even analyti- 
cally in some important cases. Although, by its design, 
the ME approach should exactly correspond to the EQ 
model when the properties of the contacts are the same, 
it is difficult to give a formal proof of the equivalence of 
the two methods because the EQ model does not have a 
known analytical solution, except in very specific cases. 
The smooth sliding, as well as the transition from stick- 
slip to smooth sliding, emerge due to a finite width of the 
distribution P c (x) and have (almost) nothing in common 
with the melting-freezing or inertia microscopic mecha- 
nisms of stick-slip. We observed that, for stick slip to oc- 
cur, the distribution of breaking thresholds P c {x) must 
have a low enough weight for small thresholds, for ex- 
ample, P c (x) oc x 3 if x — » as follows from Eqs. (jSTj) 
and (|55|) . 

The complex problem of behavior of a tribological sys- 
tem is split into two independent subproblems. In the 
present paper we studied the first one: dynamics of the 
friction contact, if the distribution of static thresholds 
P c {x) is known. The second separate problem is to find 
the distribution P c (f) for a given system. Although we 



did not focus on this problem we examined the various 
possible situations of microscopic contacts to estimate 
P c (x) for those situations in order to examine how they 
fit in the framework of the master equation approach. For 
a contact of two hard rough surfaces (the multi-contact 
interface, or MCI in notations used by Baumberger and 
Caroli [J]) this problem reduces to the statistics of as- 
perities. The contact aging, i.e. the increase of static 
thresholds with the time of stationary contact, is due to 
two processes: the first (and more important) process is 
due to geometrical aging, or the increase of contact area 
at the asperity, and the second, due to structural aging, 
or restructuring of the contact. The value of the param- 
eter D in Eq. (j4Tj) may be estimated from experiments: 
it was found that in most cases the average static thresh- 
old n s = F s / '-Fioad grows logarithmically with the waiting 
time t m , d^ s /d(\nt w ) w 10~ 2 

For flat surfaces such as, e.g., the mica surfaces used in 
the SFA experiments, the surface consists, most probably, 
of domains of different orientations (as for poly-crystals, 
and even for mono-crystals, where simply the sizes of 
domains are much larger). However, even if both sliding 
surfaces have an ideal crystalline structure on the macro- 
scopic scale, anyway the lubricant film should consist of 
domains of different orientations or different metastable 
structures separated by dislocations (domain walls). The 
parameters of these domains can be studied by standard 
methods of molecular dynamics, while their evolution 
with time (the structural aging, which should lead to 
the increase of the static thresholds) may be described, 
e.g., by the Ginzburg-Landau theory. 

It is interesting to note that the EQ-ME viewpoint dis- 
cussed here appeared in biological physics more than fifty 
year ago to describe skeletal muscles (46|. The Lacker- 
Peskin model [47l - l5fj| describes the binding-unbinding of 
the myosin molecules on actin filaments in terms of an 
equation which is very similar to our master equation. 

Finally, let us mention questions which were not con- 
sidered in the present work. First, throughout the pa- 
per we ignored inertia effects. They would imply that 
the distribution Q(x) should be extended to negative x 
values because some contacts could pick up enough ki- 
netic energy during sliding to overshoot. Returning to 
the stick (pinned) state the spring force acting on such 
a contact would be negative [14j. Moreover in the case 
of dFk(v)/dv < which appears due to contacts aging, 
the regimes of stick-slip and smooth sliding may be sep- 
arated by a regime of irregular (chaotic) motion due to 
inertia effects. Second, we completely ignored the role of 
the interaction between the contacts. The effects of in- 
teraction should change (increase) the critical velocity of 
the transition from stick-slip to smooth sliding. Interac- 
tion may be incorporated indirectly in a mean-field fash- 
ion by a renormalization of the distributions P c {x) and 
R(x). However, a concerted, or synchronized breaking 
(triggering) of contacts may be studied numerically only 
within the earthquake-like model (e.g., see Ref. Q)- It 
cannot be included rigorously in the master equation ap- 
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proach, unless it is coupled to a deformation field which 
deeply modifies the approach. Also, we considered the 
contacting planes of the sliding blocks as rigid planes, 
i.e., we ignored a possible variation (or fluctuation) of F 
in the (x, y) plane. Moreover, the latter should be cou- 
pled with the nonuniform elastic deformation of the block 
(for a preliminary work in this direction see Ref. 19]). 
Finally, we ignored a possible heating of the contacts 
due to sliding. All these question deserve separate stud- 
ies. In conclusion, note that the approach developed in 
the present paper, is to be applied to meso- or macro- 
scopic systems only and cannot be used to explain the 
AFM/FFM tip-based experiments, except if there is a 
multi-contact through a sufficient number of tip atoms. 
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line) for the rectangular distribution P c {x) (solid line). 
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One can easily find that C — x s + p 1 , 

f 1 for < x < x s 
P[X} ~ 1 e -P( x - x ^ for x > x s , 



(B2) 



Appendix A: Shear elastic constant 

To estimate possible values of the shear constant ki , let 
us assume that the contacts have the shape of a cylinder 
of radius r c and length h (i.e., h is the thickness of the in- 
terface). Also, we suppose that one end of such a contact 
column is fixed, while a force fi is applied to the free end. 
This force will lead to the displacement xi = fi/ki of the 
end, where ki = 3E c I/h 3 , E c is the Young modulus of 
the material of the contact and I = itr\ /4 is the moment 
of inertia of the cylinder (see Ref. [5l[). This leads to 

h = {3n/4){E c r c )(r c /h) 3 . (Al) 

As the Young modulus is coupled with the transverse 
sound velocity of the material by the relationship E c — 
2p(l + a)c 2 , where p is the mass density and a is the 
Poisson ratio, we obtain 

ki = (3tt/2)p(1 + n)c 2 r c {r c /hf . (A2) 

Thus, if r c « h, we obtain ki ~ pc 2 ^J~A~i, where Ai = nr 2 
is the contact area. 



Appendix B: Two simple examples of steady state 
solutions 

As a simple example, let us consider the distribution 
P(x) = p Q(x — x s ), or 

p w = {l If lit (B1) 



pn _J for < x < x s , , . 

n{X) ~ \ pe-?'*- 1 ') for x > x s , 1 6) 

and 

F k = ^N c {ki) [xs+p- 1 +p- 1 /{l+px s )] . (B4) 

In particular, if x s = 0, then Q s {x) — 0(x)pe~ px and 
F = N c (ki) /p, while for the case of x s > and p — > oo we 
obtain that Q s (x) = xj 1 within the interval < x < x s 
and outside it, so that F = ^N c (ki)x s . 

Another simple example, which admits an exact so- 
lution, is the case of the rectangular P c (x) distribution 
shown in Fig. [20l 

f if x < 0.5, 
P c (x) = < 1 if 0.5 < x < 1.5, (B5) 
[ if x > 1.5. 

The fraction density of breaking contacts P(x) is given 
byEq. ©: 

( if x < 0.5, 
P(x) = I (1.5 -x)- 1 if 0.5 < x < 1.5, (B6) 
I oo if x > 1.5. 

P{x), given by Eq. (J6j> diverges when x tends to 1.5, 
and then, according to its physical meaning it has to be 
infinite for larger values of a;. In this case 

( if x < 0.5, 

U(x) = l -ln(1.5-a;) if 0.5 < x < 1.5, (B7) 
| oo if x > 1.5, 
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E P {x) 



1.5 







if 
if 
if 



x < 0.5, 

0.5 < x < 1.5, 

x > 1.5. 



(B8) 



The earthquake model with the distribution (|B5[) and 
the master equation model with P(x) given by Eq. (|B6[) 
have exactly the same solution for any initial configura- 
tion. 

Notice that this example clearly demonstrates that 
the probability distribution P c (x), which determines the 
static thresholds for the newborn contacts, is different 
from the distribution of static thresholds P xs (x) which is 
achieved in the steady state as pointed out in Sec. [Til 



Appendix C: Friction loop 



equations for the forward and backward motions we see 
that they obey the symmetry relations P b (x) = P(—x) 
and Rb(x) = R(—x) which are a manifestation of the ir- 
reversibility of the master equation. Indeed, if the force 
fi on a given contact i approaches and overcomes f S i, 
the contact breaks; but if we now reverse the direction of 
the motion, this contact does not jump back to the value 
fi ~ fsi] instead fi decreases until it reaches the value 

fi ~ fsi- 

Although this comment does not bring in new physics, 
it is important for experimentalists. In a typical tribo- 
logical experiment, the top substrate periodically moves 
forward and backward, and as a result, the so-called "fric- 
tion loop" is observed. The same loop can easily be cal- 
culated with the ME approach as shown in Fig. [5TJ 
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Appendix D: Delay effects 



FIG. 21: Friction loop. The distribution P c (x) is Gaussian 
with x 3 = 1 and a s = 0.2, the initial distribution Qi n i(a;) is 
Gaussian with Sj n j = and <Ti n i = 0.025. The top substrate 
moves to the right, then to the left, and finally again to the 
right as indicated by arrows. 

Through the paper we assumed that the top substrate 
moves continuously to the right, i.e. AX > 0. When the 
substrate moves to the left, or AX < 0, Eqs. ([5]), ([7]) 
and ([H]) must be modified in the following manner: 



AQ_(x;X) = -P b {x)AXQ(x;X) 



(CI) 



/oo 
dx'AQ-{x'\X) (C2) 
-co 



and 



dQ(x;X) dQ(x;X) 



-P b (x)Q(x;X) 



dx dX 

= R b {x) / dx' P b (x')Q(x';X), (C3) 



where we have denoted by index b (P b (x), R b (x)) quan- 
tities relative to the backward motion. Comparing the 
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FIG. 22: (Color online): The dependences F{X) for differ- 
ent delay times AX m — (black solid), 0.2 (red dashed), 0.4 
(black dotted), 0.8 (magenta dash-dot-dotted) and 1.2 (blue 
dash-dotted curve). The distribution P c (x) is Gaussian with 
x s — 1 and cr s = 0.2, the initial distribution Qi n i(x) is Gaus- 
sian with Xi n i = 0.1 and cri n j = 0.025, and (ki) = 1. 



To establish Eq. ([7]) we assumed that, after a slip, a 
contact would stick again without any delay. Actually 
restoring the contact always requires a small delay t,„, 
particularly for the melting/freezing mechanism. Sim- 
ulation [3( shows that r TO ~ 10 2 to, where tq ~ 10~ 13 s 
is a characteristic period of atomic oscillations in the 
lubricant. Therefore to write Eq. ([7J we have to use 
AQ + (x;X + AX m ), where AX rn = vT m and v is the 
driving velocity. 

In the body of paper we ignored the delay effect. In 
this case AX m = AX, and we may drop the +AX because 
it leads to a second order correction since it appears in a 
term which is already a correction to Q. 

When the delay is taken into account, the integro- 
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FIG. 23: (Color online): The final number of asperities in 
contact (red circles) and the kinetic friction force (black di- 
amonds) in smooth sliding regime as functions of the delay 
time (other parameters as in Fig. I22[) . The solid curve shows 
the (1 + AX m ) _1 dependence. 



lution given by Eqs. (fT2TJT4|) . which describes the smooth 
sliding. However, there is one exception from this gen- 
eral scenario, when the model admits a periodic solution 
and F(X) ^ Const even in the limit X — > oo. This is 
the singular case when all contacts are identical, i.e., all 
contacts are characterized by the same static threshold 
x s , so that P c (x) — S(x — x s ), or P(x) — for x < x s 
and P(x) = oo for x > x s . 

As one can check by direct substitution, in this par- 
ticular case the steady-state solution of Eq. © has the 
following form: 



Q(x; X) = S(x - X) [G(x) - B(x - x s )} 



(El) 



where the function S(t;) is defined by the initial condition: 
S(£) = <3ini(0 f° r the interval < £ < x s , and then S(£) 
should be repeated periodically over the whole interval 

— oo < £ < oo, 



S(£) V n integer . 



differential equation © has to be modified to become In this simple case the total force © is equal to 



dQ{x;X) dQ(x;X) 



P(x)Q(x;X) 



dx dX 

r DO 

R(x) / dx' P(x')Q(x';X - AX V 



(Dl) 



The normalization condition ^ does not hold for 
Q(x; X) because the total number of unbroken contacts, 
is not conserved anymore. When an asperity breaks it is 
not in contact with the substrate during the time r m un- 
til the contact is restored. Thus, iV c varies with X, and 
N C (X) < iV c (0), because some contacts are in a tem- 
porarily broken state even at smooth sliding. Typical 
dependences F(X) for different values of the delay time 
are shown in Fig. [22l The kinetic friction force during 
smooth sliding is such that Fk oc N c (oo) and, therefore, 
it depends on pre-history of the contacts. If one starts 
from the same initial configuration, the final force Fk is 
lower when the delay time increases (see Fig. [23]) . The de- 
pendence Fk(v) is described by Eq. (|42l) with F(oo) = 

and Paging — X s /T m . 



Appendix E: The singular case 

As we have shown, in a general case the solution of the 
master equation always approaches the steady-state so- 



F(X) = N c k 



x 



(E2) 



The static friction force takes the minimal value F s = 
^N c kx s for the uniform initial distribution, Qi n i(a;) = 
xj , when F(X) does not depend on X, and the maximal 
value F s = N c kx s for the delta-function initial distribu- 
tion Qini(a;) = S(x — xq) with some < xq < x s , when 
the function F(X) has sawtooth shape changing from 
to F s . 

However, the periodic solution described above only 
exists for a distribution P c (x) with a single threshold. If 
the contacts are characterized by more than one thresh- 
old value, for example, if one part of contacts has the 
threshold x s \ and another part the threshold x S 2 ^ x s \ 
[i.e., P c {x) is described by a sum of two delta-functions], 
then the system will always approach the stationary 
steady state. This is demonstrated in Fig. [24] where we 
compare the system evolution in cases of one-peak and 
two-peak P c {x) distributions. Notice, however, that this 
statement is only valid for an infinite set of contacts (the 
number of contacts with each threshold must be infinite) 
and cannot be applied for a microscopic system where, 
e.g. two tips move over a surface. 
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